---
title: "Replication code: Who wants What, and Why? Attitude Polarization and Political Viability of UBI"
author: "Bastian Becker, Hanna Schwander"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_depth: 4
    toc_float: 
      collapsed: false
      smooth_scroll: true
    code_folding: 'hide'
    self_contained: true
---

This script replicates the results of the article 'Who want What, and Why? Attitude Polarization and Political Viability of UBI' by Bastian Becker and Hanna Schwander, forthcoming at Journal of European Social Policy. 

To ensure anonymity, the data contains open text responses only in their recoded version. 

```{r setup, include = FALSE}

# -------------------------------------------------------------
# Setup
# -------------------------------------------------------------

# Load packages
library(ggplot2)
library(tidyr)
library(dplyr)
library(cregg)

# Global settings
options(warn = -1)

# Load data
DAT <- readRDS("replication_data.Rdata")

```

# Sample descriptives

```{r descriptives}

# -------------------------------------------------------------
# Sample descriptives
# -------------------------------------------------------------

table(DAT$gender)/sum(table(DAT$gender))

table(cut(DAT$age, c(17,30,45,60,120)))/sum(table(DAT$age))

table(DAT$vote_intention)/sum(table(DAT$vote_intention))


```

# UBI Support

```{r ubi support}

# -------------------------------------------------------------
# UBI support by party
# -------------------------------------------------------------

DAT |> 
    filter(partyID_en%in%c("AfD","FDP","CDU","SPD","Green Party","The Left")) |> 
    group_by(partyID_en, ubi_support3) |> 
    summarize(n=length(partyID_en), .groups = "drop_last") |> 
    group_by(partyID_en) |> 
    mutate(n_group=sum(n))  |> 
    ungroup() |> 
    mutate(share=n/n_group*100) |> 
    ggplot(aes(y=share, fill=ubi_support3, x=partyID_en)) + 
    geom_bar(stat = 'identity', position = 'dodge') + 
    scale_fill_manual(values=c("firebrick2","darkgrey","black")) +
    xlab("") + ylab("Share by Party [%]") + 
    theme_minimal() +
    theme(legend.title = element_blank())

# -------------------------------------------------------------
# UBI support by income
# -------------------------------------------------------------

DAT$income_ind_cat<-factor(DAT$income_ind_cat, 
                           levels = c("- 750€","750-1,499€","1,500-2,499€", "2.500-3,999€", "4,000-5,999€", "6,000€ -"),
                           labels = c(c("- 750€","750 -\n1,499€","1,500 -\n2,499€", "2.500 -\n3,999€", "4,000 -\n5,999€", "6,000€ -"  ))
                           )

DAT |> 
   filter(!is.na(income_ind_cat)) |> 
    group_by(income_ind_cat,ubi_support3) |> 
    summarize(n=length(income_ind_cat), .groups = "drop_last") |> 
    group_by(income_ind_cat) |> 
    mutate(n_group=sum(n))  |> 
    ungroup() |> 
    mutate(share=n/n_group*100) |> 
    ggplot(aes(y=share, fill=ubi_support3, x=income_ind_cat)) + 
    geom_bar(stat = 'identity', position = 'dodge') + 
    scale_fill_manual(values=c("firebrick2","darkgrey","black")) +
    xlab("") + ylab("Share by Income [%]") + 
    theme_minimal() +
    theme(legend.title = element_blank())

```


# UBI Arguments

```{r arguments}

# -------------------------------------------------------------
# Argument Aggregation and Visualization
# -------------------------------------------------------------

# Define argument pairs and labels
arg_pairs <- list(
  Justice      = c("pro_justice", "con_justice"),
  Liberty      = c("pro_liberty", "con_liberty"),
  LaborMarket  = c("pro_labor", "con_labor"),
  FiscalPolicy = c("pro_fiscal", "con_fiscal"),
  WelfareState = c("pro_welfare", "con_welfare"),
  Economy      = c("pro_economy", "con_economy"),
  Democracy    = c("pro_democracy", "con_democracy")
)

# Helper function to compute argument aggregates
aggregate_args <- function(data) {
  out <- lapply(names(arg_pairs), function(arg_name) {
    pro_var <- arg_pairs[[arg_name]][1]
    con_var <- arg_pairs[[arg_name]][2]
    agg <- aggregate(ID ~ get(pro_var) + get(con_var), data, function(x) length(unique(x)))
    colnames(agg) <- c("pro", "con", "n")
    agg$arg <- arg_name
    return(agg)
  })
  
  agg <- bind_rows(out)
  agg <- pivot_wider(agg, values_from = n, names_from = c(pro, con))
  
  # Fill missing combinations
  agg[is.na(agg)] <- 0
  
  # Compute totals and shares
  agg <- agg %>%
    mutate(
      total = `1_0` + `0_1` + `1_1`,
      pro   = `1_0` / nrow(data),
      con   = `0_1` / nrow(data),
      both  = `1_1` / nrow(data)
    )
  
  return(agg)
}

# -------------------------------------------------------------
# Aggregation: All Respondents
# -------------------------------------------------------------

agg_all <- aggregate_args(DAT)

ggplot(agg_all, aes(y = reorder(arg, total))) +
  geom_bar(aes(x = pro + both / 2), stat = "identity", width = 0.8, fill = "black") +
  geom_bar(aes(x = -(con + both / 2)), stat = "identity", width = 0.8, fill = "black") +
  geom_bar(aes(x = both / 2), stat = "identity", width = 0.8, fill = "grey40") +
  geom_bar(aes(x = -both / 2), stat = "identity", width = 0.8, fill = "grey40") +
  geom_vline(xintercept = 0, color = "white", linewidth = 0.5) +
  scale_x_continuous(
    breaks = c(-0.4, -0.2, 0, 0.2, 0.4),
    labels = c(40, 20, 0, 20, 40),
    limits = c(-0.4, 0.4),
    name = "Share of All Respondents [%]"
  ) +
  ylab("") +
  labs(subtitle = "Contra                                                   Pro") +
  theme_minimal() +
  theme(
    axis.ticks = element_blank(),
    plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold")
  )

# -------------------------------------------------------------
# Aggregation: By UBI Support Group
# -------------------------------------------------------------

groups <- c("Proponents", "Opponents", "Indifferent")
agg_by_group <- lapply(groups, function(g) {
  sub <- DAT %>% filter(ubi_support3 == g)
  agg <- aggregate_args(sub)
  agg$ubi_support3 <- g
  return(agg)
}) %>%
  bind_rows()

agg_by_group$ubi_support3 <- factor(agg_by_group$ubi_support3,
                                    levels = c("Proponents", "Indifferent", "Opponents"))

agg_by_group[is.na(agg_by_group$total), c("total", "both")] <- 0

# -------------------------------------------------------------
# Plot: Aggregation by Support Group
# -------------------------------------------------------------

ggplot(agg_by_group, aes(y = reorder(arg, total), fill = ubi_support3)) +
  geom_bar(aes(x = pro + both / 2), stat = "identity", width = 0.8, position = "dodge") +
  geom_bar(aes(x = -(con + both / 2)), stat = "identity", width = 0.8, position = "dodge") +
  geom_vline(xintercept = 0, color = "white", linewidth = 0.5) +
  scale_x_continuous(
    breaks = c(-0.4, -0.2, 0, 0.2, 0.4),
    labels = c(40, 20, 0, 20, 40),
    limits = c(-0.5, 0.5),
    name = "Share per Support Group [%]"
  ) +
  ylab("") +
  labs(subtitle = "Contra                                                   Pro") +
  scale_fill_manual(values = c("black", "grey", "red"), name = "") +
  theme_minimal() +
  theme(
    axis.ticks = element_blank(),
    plot.subtitle = element_text(size = 10, hjust = 0.5, face = "bold")
  )


```


# Conjoint analysis

```{r conjoint}

# -------------------------------------------------------------
# Helper Function: Extract and Label Conjoint Sets
# -------------------------------------------------------------

extract_conjoint <- function(prefix, proposal, v16, v17_or_v18) {
  vars <- paste0(prefix, c("_1", "_2", "_3", "_4", "_5", "_6"))
  dat <- DAT[, c("ID", vars, v16, v17_or_v18)]
  colnames(dat) <- c("ID", "Citizenship", "Residence", "Jurisdiction",
                     "Requirements", "Generosity", "Financing",
                     "selected", "score")
  dat$selected <- as.integer(dat$selected == paste("Vorschlag", proposal))
  return(dat)
}

# -------------------------------------------------------------
# Data Assembly
# -------------------------------------------------------------

DAT_conjoint <- data.frame()

for (i in 1:5) {
  # A proposals
  a_set <- extract_conjoint(paste0("VGL_", i, "A"), "A",
                            paste0("v_16_", i), paste0("v_17_", i))
  DAT_conjoint <- rbind(DAT_conjoint, a_set)

  # B proposals
  b_set <- extract_conjoint(paste0("VGL_", i, "B"), "B",
                            paste0("v_16_", i), paste0("v_18_", i))
  DAT_conjoint <- rbind(DAT_conjoint, b_set)
}

# -------------------------------------------------------------
# Variable Recoding
# -------------------------------------------------------------

DAT_conjoint$score <- as.numeric(factor(
  DAT_conjoint$score,
  levels = c("(0) Auf keinen Fall", "(1)", "(2)", "(3)", "(4)",
             "(5)", "(6)", "(7)", "(8)", "(9)", "(10) Voll und ganz")
)) - 1

DAT_conjoint <- DAT_conjoint %>%
  mutate(
    Requirements = factor(Requirements,
      levels = c("Keine", "Einkommen unter Existenzminimum",
                 "Vorübergehende Arbeitslosigkeit",
                 "Gemeinschaftsarbeit oder Pflegetätigkeit"),
      labels = c("None", "Low income", "Temporary unemployment", "Community/care work")
    ),
    Citizenship = factor(Citizenship,
      levels = c("Alle", "EU-Staatsbürger", "Deutsche Staatsbürger"),
      labels = c("All", "European Union", "German")
    ),
    Residence = factor(Residence,
      levels = c("Keine", "Ein Jahr", "Fünf Jahre"),
      labels = c("0 years", "1 year", "5 years")
    ),
    Jurisdiction = factor(Jurisdiction,
      levels = c("Deutschland", "Hamburg"),
      labels = c("Federal", "State")
    ),
    Generosity = factor(Generosity,
      levels = c("550€ (angelehnt an Bürgergeld-Regelsatz)",
                 "1100€ (angelehnt an Existenzminimum)",
                 "1700€ (angelehnt an Mindestlohn bei Vollzeitbeschäftigung)",
                 "2300€ (angelehnt an Durchschnittseinkommen)"),
      labels = c("550EUR", "1,100EUR", "1,700EUR", "2,300EUR")
    ),
    Financing = factor(Financing,
      levels = c("Kürzung anderer Sozialleistungen", "Aufnahme von Schulden",
                 "Erhöhung der Einkommenssteuer", "Einführung einer Vermögenssteuer"),
      labels = c("Welfare cuts", "Public debt", "Income tax", "Wealth tax")
    )
  )

# Add UBI support
DAT_conjoint <- merge(DAT_conjoint, DAT[, c("ID", "ubi_support3")], by = "ID")
DAT_conjoint$ubi_support3 <- as.factor(DAT_conjoint$ubi_support3)

# -------------------------------------------------------------
# Analysis: Marginal Means (Choice)
# -------------------------------------------------------------

mod_main <- cj(DAT_conjoint,
  selected ~ Generosity + Requirements + Citizenship + Residence + Financing + Jurisdiction,
  id = ~ID, estimate = "mm"
) %>% mutate(across(c(estimate, lower, upper), ~ .x * 100))

mod_cond <- cj(DAT_conjoint,
  selected ~ Generosity + Requirements + Citizenship + Residence + Financing + Jurisdiction,
  id = ~ID, by = ~ubi_support3, estimate = "mm"
) %>% mutate(across(c(estimate, lower, upper), ~ .x * 100))

mod_main$mod <- "Average Support"
mod_cond$mod <- "Conditional Support"
mod_main$ubi_support3 <- "All"

common_cols <- intersect(colnames(mod_main), colnames(mod_cond))
mod_combined <- rbind(mod_cond[, common_cols], mod_main[, common_cols])
mod_combined$ubi_support3 <- relevel(factor(mod_combined$ubi_support3), ref = "All")

# Plot Marginal Means
plot(mod_combined, feature_headers = FALSE, group = "ubi_support3",
     vline = 50, vline_color = "darkgrey") +
  aes(shape = ubi_support3, fill = ubi_support3) +
  scale_shape_manual(values = c(4, 25, 18, 24)) +
  scale_colour_manual(values = c("black", "red", "grey", "black")) +
  scale_fill_manual(values = c("black", "red", "grey", "black")) +
  coord_cartesian(xlim = c(37, 63)) +
  xlab("Selection probability [%]") +
  theme_minimal() +
  theme(legend.position = "bottom") +
  facet_grid(feature ~ mod, scales = "free_y")


# -------------------------------------------------------------
# Analysis: AMCE (Choice)
# -------------------------------------------------------------

mod_main_amce <- cj(DAT_conjoint,
  selected ~ Generosity + Requirements + Citizenship + Residence + Financing + Jurisdiction,
  id = ~ID, estimate = "amce"
)

mod_cond_amce <- cj(DAT_conjoint,
  selected ~ Generosity + Requirements + Citizenship + Residence + Financing + Jurisdiction,
  id = ~ID, by = ~ubi_support3, estimate = "amce"
)

mod_main_amce$mod <- "Average Support"
mod_cond_amce$mod <- "Conditional Support"
mod_main_amce$ubi_support3 <- "All"

common_cols <- intersect(colnames(mod_main_amce), colnames(mod_cond_amce))
mod_combined_amce <- rbind(mod_cond_amce[, common_cols], mod_main_amce[, common_cols])
mod_combined_amce$ubi_support3 <- relevel(factor(mod_combined_amce$ubi_support3), ref = "All")

# Plot AMCE
plot(mod_combined_amce, feature_headers = FALSE, group = "ubi_support3",
     vline_color = "darkgrey") +
  aes(shape = ubi_support3, fill = ubi_support3) +
  scale_shape_manual(values = c(4, 25, 18, 24)) +
  scale_colour_manual(values = c("black", "red", "grey", "black")) +
  scale_fill_manual(values = c("black", "red", "grey", "black")) +
  theme_minimal() +
  theme(legend.position = "bottom") +
  facet_grid(feature ~ mod, scales = "free_y") +
  xlab("Estimated AMCE")


# -------------------------------------------------------------
# Analysis: Marginal Means (Score)
# -------------------------------------------------------------

mod_main_score <- cj(DAT_conjoint,
  score ~ Generosity + Requirements + Citizenship + Residence + Financing + Jurisdiction,
  id = ~ID, estimate = "mm"
)

mod_cond_score <- cj(DAT_conjoint,
  score ~ Generosity + Requirements + Citizenship + Residence + Financing + Jurisdiction,
  id = ~ID, by = ~ubi_support3, estimate = "mm"
)

mod_main_score$mod <- "Average Support"
mod_cond_score$mod <- "Conditional Support"
mod_main_score$ubi_support3 <- "All"

common_cols <- intersect(colnames(mod_main_score), colnames(mod_cond_score))
mod_combined_score <- rbind(mod_cond_score[, common_cols], mod_main_score[, common_cols])
mod_combined_score$ubi_support3 <- relevel(factor(mod_combined_score$ubi_support3), ref = "All")

# Plot Score-Based Marginal Means
plot(mod_combined_score, feature_headers = FALSE, group = "ubi_support3",
     vline_color = "darkgrey") +
  aes(shape = ubi_support3, fill = ubi_support3) +
  scale_shape_manual(values = c(4, 25, 18, 24)) +
  scale_colour_manual(values = c("black", "red", "grey", "black")) +
  scale_fill_manual(values = c("black", "red", "grey", "black")) +
  theme_minimal() +
  theme(legend.position = "bottom") +
  facet_grid(feature ~ mod, scales = "free") +
  xlab("Profile score [0–10]")

```

